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Chapter 0 


INTRODUCTION 


OBJECTIVE: 

The objective of this project was to enhance the computational fluid dynamics 
(CFD) code NEKTON to include additional capabilities necessary for chemical vapor 
deposition (CVD) applications. The major task associated with this was to mod- 
ify the incompressible code to include compressibility. This was important because 
compressibility is essential to the accurate solution of CVD problems. Other tasks 
related to CVD, and specifically, microgravity CVD, included addition of arbitrary 
body forces for g-jitter, addition of variable properties for temperature-dependent and 
species-dependent materials properties, and the ability to calculate chemical reactions 
in multicomponent reactions. 

IMPLEMENTATION : 

The strategy for the implementation of the capabilities necessary for CVD ap- 
plications was to introduce sequentially each feature necessary for solution of CVD 
problems. Each feature added to the code was tested on sample problems to verify 
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its accuracy. 

Nektonics is currently planning to release a NEKTON/CVD code in the context 
of its next major product release (Nekton 3.0). 

ORGANIZATION: 

In the following we detail chronologically the implementation of each individual 
feature in the form of summaries of the work done in each quarter. 

In Chapter 1 we describe the implementation and demonstrate the use of the 
capability for solving problems involving the “g-jitter” forces. These time-varying 
forces typically are large in relation to the constant accelerations in microgravity 
environments. 

In Chapter 2 we demonstrate the use of the variable diffusivity of chemical species 
in NEKTON. We demonstrate the effects of variable diffusivity on mass transfer 
problems in geometries prototypical in CVD applications. 

In Chapter 3 we describe the multibody radiation capability in NEKTON. We 
describe the module which gives the additional capability of calculation of Gebhart 
factors through inversion of the shape factor matrices. 

In Chapter 4 we show the implementation of the Soret effects in the passive scalars 
of NEKTON. 

In Chapter 5 we describe the implementation of the Motif graphical user interface 
in NEKTON. We give examples of the use of it for ease of the input of geometries 
and ease of evaluation of the results. 

In Chapter 6 we first give the expansion of the Navier-Stokes equations using 
perturbations for small Mach numbers. We justify the use of the resulting Low-Mach 
number approximation for CVD applications and demonstrate its use in compressible 
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test cases. 


In Chapter 7 we present the method of Bussing and Murman for the solution of 
coupled problems in which the time scales for chemical reactions was much smaller 
than the diffusive time scales. We demonstrate the use of this point-implicit method 
in test problems with varying time scales. 

In Chapter 8 we show the method for solving the convection-diffusion problem for 
nondilute mixtures of gases. This involves calculation of a convective component of 
the diffusion process and a mixture-dependent density. 

Finally, in Chapter 9 we summarize the work done in this project. We discuss the 
physical processes important in CVD problems, the approach we used to model these 
processes, and the results we have achieved. 
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Chapter 1 


Unsteady Body Forces (g-jitter) 


OBJECTIVE: 

The technical objective to be attained by the close of the first quarter in Phase II 
was to develop the capability to handle arbitrary fluctuating three-dimensional body 
forces (”g-jitter”). 

BACKGROUND: 

The motivation for enhancing the simulation capabilities of NEKTON through 
inclusion of an arbitrary three-dimensional body force stems from the need to ac- 
curately handle the prevailing gravity levels in space during materials processing. 
The recorded ^-levels on board the shuttle appear to establish that the direction and 
magnitude of the gravity vector in space vary with time [1]. Whereas fundamental 
considerations suggest a minimum residual gravitational acceleration of 10 -6 of earth 
gravity in space, current measurement techinques appear to be limited to 10 -5 g e . 
Further more, the lowest frequency of the ^-jitter is found to be close to 20 Hz. As 
the primary driving force for interference of convection with crystal growth in space 
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is associated with the buoyancy forces in the reactor, it is reasonable to expect that 
convection in these systems is a complex function of time leading to three-dimensional 
flow structures. Numerical simulation of crystal growth processes in space (e.g. [2]) 
are generally based on a uniform magnitude and direction of g. The time-dependant 
aspect of gravity level has been considered by several researchers [3-6]. It has been 
shown in ref. [6] that the response of convective systems to time-dependant gravity 
is controlled by the ratio of vibrating frequency of g to the diffusive time scale in 
the system, the so-called Wormerley number W = h 2 u>/v. For W < 1 convection in 
the system follows the time-dependant change in g, and for W > 1 the intensity of 
convection in the melt decreases with W. The formulations developed in refs. [3-6], 
amongst many others, represent idealizations of the prevailing gravity field in the 
space. Recently Alexander and colleagues have shown that impulses in g as well as 
the three dimensional nature of the gravity field can significantly influence the dopant 
distribution during directional solidification of doped semiconductors [7]. With the 
continuing advances in accelerometer technology reliable data on ^-levels in space can 
be expected to be available in the near future. As the magnitude and direction of 
g determines the flow structure in crystal growth systems, the presently developed 
capability of NEKTON can provide for accurate simulation of convection in crystal 
growth systems using measured unsteady 3-dimensional (/-levels in space. 

IMPLEMENTATION . 

The gravitational body force was input in the right hand side of the Navier stokes 
equations in a manner analogous to the handling of the nonlinear convective term. 
That is, a third order Adams-Bashforth integration scheme was used. This scheme, 
in addition to providing third order accuracy for the integration of a time-dependant 
gravitational body force, has favorable stability characteristics. No additional stabil- 
ity constraints beyond the Courant condition were imposed by the oscillating body 
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force. The user accesses this gravitational body force via 3 FORTRAN functions 
typed in during problem definition in the preprocessor. Through each of the 3 FOR- 
TRAN functions the user can define 3 independent components of the gravitational 
force, each of which can be either constant or can vary with time according to an 
arbitrary functional relationship. 

The body force was tested using three cases. In each case the steady component 
of the Rayleigh number, that is, the magnitude of the convective force caused by 
the steady component of gravity, was 10. In cases b and c a sinusoidally oscillating 
component of component of gravity, resulted in the oscillating Rayleigh number of 
100. The problem was solved for the hemispherical geometry depicted in figure 1. A 
temperature of 100 was imposed at the seed surface, and a temperature of zero was 
imposed at the inner surface of the support shaft. 


Case a 

Wormersley 

Number 0 

Rayleigh 

Number 

Steady 10 

Oscillating 100 


Case b Case c 


0.1 100 


10 

10 

100 

100 


Here we define the Wormersley number as W=h**2 omega/nu, where h is the 
sphere radius, omega is the frequency in cycles/time, and nu is the kinematic viscos- 
ity. The test cases represent a zero Wormersley number control, a low Wormersley 
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number, and a high Wormersley number. 

For relatively low Wormersley numbers (case b), we expect the dynamic compo- 
nents of the oscilliating force to be of marginal unimportance; we expect a quasi- 
steady flow at each instant corresponding to the steady flow resulting from a steady 
gravitational force equal to the gravitational force at that instant. 

In the high w limit (case c), we expect a strong interaction between the dynamic 
component of the oscillating force and the dynamics of the flow regime. In this case, 
w’e expect a filtering behavior such that the response of the system in terms of velocity 
magnitude to decrease with Wormersley number (frequency). 

Note that these test cases were run each with a single frequency component of 
the gravitational force. This in no way reflects a code limitation; the structure of 
the function input allows for multiple frequencies and directions to be simultaneuosly 
input. Studies of nonlinear interactions between frequency components are therefore 
possible with this implementation. 

RESULTS: 

In figures a, b, and c we plot the axial velocity component at a point on the axis 
about midway between the crystal surface and the hemisphere surface. In case a the 
velocity reaches a peak of 0.33e-2 before reaching its steady value of approximately 
O.le-2. This velocity overshoot is caused by the evolution of the temperature field 
from its initial conditions. As steady state is reached the convection set up by the 
velocity reduces the temperature difference within the fluid, acting to relieve part 
of the force driving the fluid. In case b the velocity reaches a steady periodic half- 
amplitude of O.le-2. In case c the velocity reaches a steady periodic half-amplitude of 
0.12e-3. There is a clear indication of suppression of the effects of convective force with 
increasing Wormersley number. For Wormersley number of 0.1, the magnitude of the 
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oscillating gravitational component needed to be an order of magintude higher than 
the steady case to achieve com- parable convective velocity magintudes. Increasing 
the Wormersley number to 100 supresses the convective velocity by an additional 
order of magnitude. 

REFERENCES 

1. N. Trappen and F.J. Demond, Porceedings of the Nordemey Symposium on 
Scientific Results of the German Spacelabd Mission Dl, Norderney, Germany, 
1986. 

2. P.M. Adornato and R.A. Brown, J. Crystal Growth, 80 (1987) 155 

3. LAV. Spradely, SAV. Bourgeois, and F.N. Lin, AIAA paper, 75-695, 1975. 

4. Y. Kamotani, A. Parsad, and S. Ostrach, AIAA Journal, 19 (1981) 511. 

5. P.R. Griffin and S. Motakef, Applied Microgravity Technology, II (1989) 121- 
132. 

6. J.I.D Alexander, J. Ouzzani and F. Rosenberger, J. Crystal Growth, 97 (1989) 
285. 
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Chapter 2 


Variable Diffusivity 

OBJECTIVE: 

The objective in quarter 2 was to complete implementation of variable diffusivity 
of the chemical species within the context of an incompressible NEKTON. We have 
completed the implementation of this capability in NEKTON. We have run two sim- 
ulations that demonstrate the effects of variable diffusivity in a prototypical CVD 
reactor. These simulations also serve as a preliminary invesigation on the effects of 
variable properties in prototypical CVD applications. 

OTHER PROGRESS: 

Parallel to the development of new capabilities within the NEKTON solver is the 
development of the user interface. During this quarter we addressed two aspects of 
user interface development. 

First, we expanded the capability of the PRENEK preprocessor to enable the user to 
access the new features added to the NEKTON solver. PRENEK now prompts the 
user for variable or constant properties within the passive scalar fields. This enables a 
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user to enter, from the keyboard, a FORTRAN function which defines the diffusivity 
for the chemical species stored within the passive scalar fields. This allows the species 
diffusivity to vary arbitrarily with temperature. 

Second, we enhanced the structure and layout of the user interface to make the 
package a more useful CVD tool. A modern CVD tool should have an interface which 
is flexible, and readily accessible to its user without requiring extensive training or 
programming skills. Toward this goal we have adapted PRENEK and POSTNEK 
to run through the Motif toolkit for X Windows applications. In figure 1, we show 
the layout and structure of the Motif- based POSTNEK application. In this figure 
two windows are brought by POSTNEK, a window for flow control and text display 
at the top of the figure, and the plotting window below. The top level control of 
POSTNEK is handled by a menu bar which spans the top of the control and text 
wundow. Within the menu bar are sets of cascading pull-down menus which give a 
top-level, intuitive, and standardized system for control of POSTNEK. 

As many MOCVD applications can be treated as quasi-steady problems, we have 
implemented the steady state solution technique of Sidi and Celestinafl] in NEKTON. 
The advantages of this addition are twofold. First, solution of the problem requires 
less input from the user to define the problem. Second, significant savings in cpu 
time can result from accelerating the convergence to steady state. At the expense of 
losing accuracy during the transient phase, the accurate steady state solution can be 
obtained much more quickly. In appendix A we present a more detailed description 
of this technique. 
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VARIABLE DIFFUSIVITY: 


BACKGROUND: Constant gas properties are assumed within CVD applications 
in a first examination in order to make the problem more amenable to theoretical 
analysis and to make numerical investigation simpler. Real CVD applications can 
have gases in which absolute temperatures vary by a factor of two or more, which 
result in diffusivities changing by a factor of 3-4. 

IMPLEMENTATION : 


We have implemented in NEKTON the ability to input gas properties (diffusivity 
within passive scalars) which can vary locally with temperature and passive scalar 
concentration. In NEKTON this has required us to change the equation solved for the 
diffusion operator for passive scalars. Whereas we previously calculated the term as 
(D*laplacian (concentration), we now calculate the gradient(D*gradient(concentration)). 
In the latter calculation, the diffusivity D can vary arbitrarily with space without 
compromising the accuracy of the result. 

TEST CASES 


We have run two test cases. The purpose of these tests was to demonstrate the 
impact of variable properties on quantities of interest in a realistic CVD geometry, 
i.e., concentrations profiles of chemical species and deposition rates. 


In test case a we simulate diffusion of TMGa in H 2 using a constant diffusivity based 
on the bulk temperature of the ambient gas. 


In test case b we simulate diffusion of TMGa in H 2 using the dependencies from 
Moffat and Jensen [2]. This diffusivity varied according to the relation 


D = 


2.23 x 10 -5 T 173 


P 


cm 2 /sec 
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where D is in cm**2/sec and p is in atm. In figure 2 we show the geometry of the 
hemispherical reactor in which we perform our simulations. In this geometry the 
seed material temperature is fixed at 1300 K, while the growing crystal is cooled 
from below and has a somewhat lower surface temperature. In figure 3 we plot this 
temperature along line segment a-b defined in figure 4, which runs perpendicularly 
from the substrate surface toward the reactor wall. In figures 5a and 5b we plot 
along the same line segment the concentration profiles for test cases (a) and (b), 
respectively. The profiles for the two cases are qualitatively similar. As expected, 
the concentration gradient (and therefore, deposition rate) at the surface is steeper 
in case (b). This follows from the lower local diffusion coefficient in case (b) due to 
the lower than average local temperature. The net difference in deposition rate is 
approximately 25 percent. 

REFERENCES 

1. Moffat and Jensen, Journal of Crystal Growth, 77 (1986) 108-119. 

2. Sidi, Avram, and Celestina, Mark “Convergence Acceleration for Vector Se- 
quences and Applications to Computational Fluid Dynamics”, NASA Technical 
Memorandum 101327, ICOMP-88-17. 
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APPENDIX: STEADY STATE ACCELERATION 
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Steady State Acceleration 

NEKTON uses a semi-implicit approach to solve the unsteady 
Navier-Stokes equations (or unsteady convection-diffusion 
equations). If the solution has a steady state, this solution is 
obtained by a time-integration of the unsteady equations until 
steady state is reached. 

For larger Re-numbers (or Pe-numbers), the number of times 
steps, N t = x/A t , can be large. The reasons for this are: (1) the typical 
time (x) to reach steady state is proportional to the Re-number; (2) 
larger Re-number problems generally require more resolution 
(smaller mesh size, A*); (3) the time step restriction due to the 
Courant criterion, A t < CU/A* , is more severe for higher Re-numbers 
and high resolution meshes. Here U is a characteristic velocity. 

In order to speed up the steady Navier-Stokes solver, we have 
recently implemented a steady state acceleration procedure based on 
the work of Sidi et al ["Convergence Acceleration for Vector 
Sequences and Applications to Computational Fluid Dynamics", 

August 1988]. The main idea is to compute the initial transient 
behavior, and then apply an extrapolation procedure based on a 
sequence of (global) solution vectors. Figure 6 shows an example 
where this strategy was applied, resulting in a speed-up of a factor 
of 5. 
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Chapter 3 


Multibody Radiation 


OBJECTIVE: 

The objective for this report period was to implement multibody radiation capabil- 
ity in NEKTON and demonstrate its effects on the thermal calculations in a CVD 
rotating pedestal reactor geometry. 

BACKGROUND: 

The high temperatures found inside CVD reactors ensure that radiation will play a 
dominant role in heat transfer. Therefore, radiation heat transfer cannot be ignored in 
CVD reactors. Typically, radiation is handled in numerical code with the assumption 
of single-body radiation, that is, independent radiation exchange between a single 
isolated body and a nonreflecting environment. This approach breaks down when 
radiative interaction between components is important, when the environment tem- 
perature is unknown, and when the environment cannot be treated as an isothermal 

black body. 

Calculation of multibody radiative exchange allows for the calculation of tempera- 


29 



tures and heat exchange rates between multiple bodys with unknown temperatures, 
arbitrary emissivities and arbitrary relative configurations (shape factors). 

IMPLEMENTATION : 

As far as we have been able to determine, multibody radiation has not yet been 
implemented in spectral code. To implement this feature in NEKTON, the goal was 
to retain the spectral accuracy inherent in the rest of the calculations, to maintain a 
reasonable operation count, to implement it in a manner that is clean and amenable 
to easy future maintenance of the code, and to make its use relatively easy for the 
user. 

The spectral accuracy was retained by using an integration technique for the radiative 
heat flux calculation consistent with the spectrally accurate integration used in other 
parts of the code. The integrals of heat transfer between elemental edges were done 
as follows. The temperatures used were those at each Gauss-Lobatto collocation 
point on each elemental edge involved in the multibody radiation calculation. The 
area associated with each collocation point for the calculation of shape factors was 
weighted by the mass matrix. Posing the problem in this manner enabled coupling 
the spectral code with standard methods for the calculation of radiation heat transfer 
(shape factors, etc.). The resultant operation count, therefore, was the same (per 
point) as that for radiation implemented in a lower order finite difference or finite 
element method. 

The shape factors for each point was calculated based on the effective area associated 
with each collocation point. Each of these areas was set to be a panel and (in 2-D) 
Hottel’s Methods is used to calculate the shape factors. For the axisymmetric cases, 
the values from a set of standard axisymmetric configurations were adapted to the 
general case using the algebra governing shape factors. Thus, shape factors are cal- 
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culated automatically for 2-d or axisymmetric simple (without blocking) enclosures. 
For more complex systems with blocking of radiation or three dimensional problems, 
a subroutine is provided in which the shape factors can be input directly by the user. 

We chose the enclosure analysis method of Gebhart [1] for the radiative coupling 
calculation. In this approach, the local emissivities and shape factors are converted 
into Gebhart factors, which are essentially Green's functions for calculating the net 
radiative flux to each node due to the sum of the other radiative fluxes. 

Calculation of each column [G*] in the matrix of Gebhart factors is as follows: 


G lk = G-k‘k + G-iP> c u + G- 2 P 2 G 3 * + 

4- G-kPkG\k + ■ ’ • + G-.vP.vC.vic 

G lk = G - k*k + G-iPi^u +■ + 

+ G-kP/'k* + • • ■ + G -NpsGsk 


Gsk =* G-k«k + G-lPl C lk + G-2P2 G 2k + 

+ G-»PkGkk- + ■ • ' + G -hPhG»„ 


where Fj-k * s the shape factor between points j and k, pk and e*, are the reflectivity 
and emissivity, respectively, at point k. For column k of [G], the relation is expressed 
in matrix notation as: 


[G k ] = e fe [F fc ] + [F}diag[p}[G k ] 


Defining [H] = [F] diag[p] 
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(1) becomes 

-**[**] = [H] [G k ] - [G k ] 


or 

-£<,[Fjbl = ([H].[I]) [GJ 

so that 

[G t ] = ([ff| - [/])-' WFi]) 

([//’] — [I])~ l is inverted using standard direct (LU) solvers. The direct method is more 
appropriate here as the matrix is small (it contains only the radiative surface points) 
and is full. Also, the inverse needs to be calculated only once before performing N 
back solves to calculate N Gebhart vectors. 

Once the Gebhart factors are calculated, we use these to calculate the heat flux 
between the radiating surfaces: 

Q k = A„€ k aT( - ZU A,e k <,G, k T? 

where the first term on the RHS represents heat radiated from surface k and the 
second term represents the heat radiated to surface k from the N surfaces. Here a is 
Boltzmann’s constant and A k is the area of surface element k. 

This coupled set of nonlinear equations is solved iteratively in conjunction with the 
conduction heat transfer. The first term is linearized and handled semi-implicitly, 
similar to the manner in which single body radiation is calculated: 

A k t k crTk = {A k e k <7Tk)Tk 

where the h T = (Aktk&T%) has the form of a nonlinear convection coefficient which is 
updated during the iteration procedure. 
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The second term needs to be handled explicitly as a source term on the right hand 
side of the heat equation, so the procedure becomes, looping over m iterations: 

[V 2 + /C]7< m+1) = F r £f =1 Aje k aG jk {T™y 

While the implicit integration of the first term guarantees unconditional stability, the 
explicit handling of the sum on the right hand side of the equation does not. In fact 
there is a relaxation factor 0 < F r < 1 applied in the second step 

^n(m + 1) rpm p 

which, unfortunately, is problem-dependent and is adjusted to govern the stability 
and convergence rate. 

TEST CASES 

To demonstrate the multibody radiation capability and to show its effects on a CVD 
reactor, we ran three cases with the geometry of the rotating pedestal reactor shown 
in figure 1. These demonstrate the competing effects of multibody radiation (figure 
2) vs. conduction in the reactor. In each case we fix the temperature at input and 
input a fixed heat flux just below the pedestal surface. We nondimensionalize using 
the reactor diameter as the length scale, the inlet temperature as the nondimensional 
temperature, and the heat flux originating at the pedestal surface. 

The radiation cavity consists of the disk shaped surface at flow inlet, the cylindrical 
reactor wall and the disk shaped pedestal surface. These are modelled using the 
elemental edges II and 12, W1 and W2, and P, respectively (see figure 3). Each of 
these 5 elemental edges has, in turn, 5 internal collocation points associated with its 
fourth-order polynomial. Each collocation point has an effective area based on its 
spectral weighting function. 
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A typical velocity field through such a cavity is plotted in figure 4. 

The calculation of shape factors is done within the multibody radiation module of 
NEKTON based on the coordinates of each collocation point and its normal vector. 
This 25x25 shape factor matrix is inverted according to equation [1] to generate the 
Gebhart factor matrix. This matrix is full, including diagonal entries, since with 
nonunity emissivities each node can “see” the effects of its own emittance through 
reflections from other surfaces. Thus, for this problem, 625 radiation exchange paths 
are followed. 

Note that the radiation boundary conditions are superimposed upon the other bound- 
ary conditions; this allows flexibility in combining radiation with other mechanisms 
of heat transfer. The heated pedestal, for example, inputs a net of one unit of flux 
irrespective of the radiational and conduction exchange within the cavity. 

We run three cases with this geometry and set of boundary conditions. In case emissO, 
the emissivities of all radiating surfaces is set to 0. In case emiss05, the emissivities 
are set to 0.5. In case emissl the emissivities are set to 1. 

The case with emissivity of 0.0 (figures 5a, 5b, 5c) is effectively a conduction problem, 
as no radiational exchange can occur. 

The case with emissivity of 0.5 (figures 6a, 6b, 6c) exercises the full multibody radi- 
ational capability with multiple reflections between all surfaces. 

The case with emissivity of 1.0 (figures 7a, 7b, 7c) is effectively a set of single body 
radiational exchanges, as no reflections are allowed to occur. 

The effects of the multibody radiation on the resultant temperature on the overall 
solution can be measured by the the peak temperature. This peak temperature ap- 
pears near the center of the heated pedestal. It is reduced from 1.58 to 1.26 to 1.19 
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as the emissivity is increased from 0.0 to 0.5 to 1.0, demonstrating the effects of the 
increasing influence of radiative heat transfer. Moreover, the increased radiational 
component of the heat transfer causes a cooling of the center for more uniform tem- 
perature profile across the pedestal surface, as the center of the pedestal has the most 
direct view of the coolest radiating surface, the inlet. 

Note that in figure 5c, an essentially linear wall tempreature profile reflects the nature 
of the solution to the conduction problem. In figure 6c and 7c we see the devitation 
from this as the effects of the fourth power dependence of radiation and the direct 
interaction between the distant radiative surfaces causes a steeper increase in tem- 
perature profiles immediately downstream of the inlet. 

REFERENCES 

1. Gebhart, B: Surface Temperature Calculations in Radiant Surroundings of Ar- 
bitrary Complexity- For Gray, Diffuse Radiation, Int. J. Heat Mass Transfer, 
vol. 3, no. 4, pp 341-346, 1961. 

2. Seigel and Howell, Thermal Radiation Heat Transfer, Appendix E. 

3. Holman, J.P. Heat Transfer, McGraw-Hill, 1976. 
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Chapter 4 


Soret Effects 


OBJECTIVE: 

The result achieved in this report period was to implement the Soret effect capability 
for passive scalars in NEKTON and to demonstrate its effects on the chemical mass 
concentration species for calculations in a CVD rotating pedestal reactor geometry. 

BACKGROUND: 

The chemistry in the gas phase in CVD reactors and, therefore, the concentration of 
each species of mass involved in the reaction is important in finding the overall rates of 
deposition and the spatial uniformity of that deposition. The high temperatures and 
large temperature gradients inside CVD reactors cause a coupling effect between the 
heat transfer and the mass diffusion to become significant. The effect of temperature 
gradients on mass flux is known as “thermal diffusion” or the “Soret Effect”. The 
inverse effect of mass concentration gradients on heat flux is known as the Dufour 
Effect”. In MOCVD calculations, [1,2] the Soret effect is significant in effecting the 
deposition rate while the Dufour effect is not. We therefore concentrate our efforts 
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on implementing the Soret effects in order to increase the accuracy of our MOCVD 
calculations. 

IMPLEMENTATION : 

The implementation in NEKTON is valid for two- and three-dimensional as well 
as axisymmetric calculations. As most MOCVD applications are in axisymmetric 
geometries, we will use the axisymmetric form of the convection-diffusion equations 
to illustrate the implementation of these in NEKTON. 

The mass balance and associated boundary conditions for the transport of chemical 
species are, then: 


^ TMG 


i 

3 *tmg 

^TMG 



3/ 

+ < 

r 

3r 

+ B ' 3; ) 



_ 1 

3 j 

f 

cD 

\. 


TMG r 

9 In 

T 

r 

3 r 1 

tmg.h/ 

3, 

dr 



d_ 



TMG , 

3 In 

T 

+ 

d: 

DMG H, 

3: 

9r 



(7) 


'TMG.H. 


^tmg , , 9 In r 

9; kj ~Tr 


kx 


TMG 


surface 


( 8 ) 


Here D is diffusivity, c is concentration, x is mass fraction, k is a constant. The 
subscripts TMG, and H 2 represent the species to which the subscripted variable refers. 

The Soret effect consists of the term of the form 

~ 7 • ( K ^ tp) = f ik) _TL 
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where \{ is a constant, 'f is a scalar field, and f is a forcing term, here 
the Laplacian of the log of the temperature. Here this is expressed in the “strong” 
form. Using the standard finite element approach, we convert this to the weak form 
by multiplying by a test function , and integrating over the domain JI_ ■ 

^ _ y . (j<v <p)\r ^ Si. ’ ^ f M -a_ 

_A_ 

Integrating the left hand side by parts, yields 

- t ^ TT <(s V £ k S? f v v Asi. ’ 

-xu ^ 

Note that the surface (boundary) term on ^ JL fro™ this integration is exactly 
equal to the additional boundary condition imposed by the original strong statement 
of the Soret equation. Thus, imposition of the variational approach and its “natural” 
boundary conditions handles the Soret term very consistently. 

TEST CASES 

We demonstrate Soret effect with two test cases. The first test case is in a very simple 
geometry. This is to demonstrate the mechanics of the operation of the program and 
to have a result simple enough to be amenable to theoretical analysis. 

In case I, we impose a temperature onto the concentration such that the log(T) is the 
parabolic function log(T) = (X-X**2)/2 of figure la. The Soret effects constant is set 
at 1.0. The initial concentration was C=1.0. The boundary conditions were no mass 
flux on all boundaries. The simulation was run out to a nondimensional time of 10, 
which was lOx the characteristic diffusion time of the system. 
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These conditions resulted in the steady concentration plotted in figure lb. An integra- 
tion of the concentration yields a total mass in the domain within 1% of the original 
mass, with the discrepancy attributable to time stepping errors. While the concen- 
tration gradients are nonzero at the boundaries, implementation of the associated 
boundary conditions yields no net flux. This is consistent with the no-flux boundary 
conditions imposed and allows the concentration to in fact reach this steady state. 
This is required for simulations to be physically realistic. 

The quantitative comparison requires that for steady state, the dC/dx = - d(log(t))/dx. 
Comparision of the plots la and lb indicate that C and log(t) are parabolic profiles 
whose magnitude are both 0.125 and whose orientation are reversed, which verifies 
this requirement. 

The second test case is in the axisymmetric rotating pedestal CVD reactor geometry 
of figure II. Here we impose the dimensionless temperature of 1.0 at inlet and flux 
of 1.0 at the pedestal surface. The resultant temperature profile we plot in figure III 
along A-B is essentially linear. We impose the value of 1.0 concentration at the inlet 
and 0.0 at the pedestal surface. In figure IVa we plot the concentration contours. In 
figure IVb we plot the profile of concentration along segment A-B. Note the essentially 
linear nature of the concentration profile where diffusion is the primary mechanism 
of mass transfer. 

In figures Va and Vb we plot the results of the same calculation with the Soret 
effects enabled. The thermal diffusion is driven by the temperature field calculated in 
the heat transfer module of NEKTON. Here the essentially linear temperature profile 
results in log(T) convex in x, giving a concentration profile with the expected concave 
profile in x of figure Vb. 
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Chapter 5 


User Interface 


OBJECTIVE: 

The result achieved in this report period was to consolidate the implementation of 
the enhancements to the NEKTON user interface in order to accomodate the features 
(variable properties, need for easy axisymmetric and 3-D geometrical input, fast color 
fills, true 3-D scrollbar rotation, and scrollbar zoom) needed for CVD calculations and 
to provide enhanced ability to evaluate the results. 

BACKGROUND: 

Modern CAD and CAE packages rely increasingly on Graphical User Interfaces 
(GUI’s) to support increasingly complex tasks done by the user. The object of these 
GUI’s is to help the software developer produce a package that is intuitive, powerful, 
and easy to use. These GUI’s also help to standardize the interfaces by providing 
guidelines for developers that result in coordination in the development of widely 
varying packages [1-7]. Thus, a user familiar with editting a file on a Macintosh or 
with DecWindows will recognize some of the essential features in panning, zooming, 
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and pulldown menu selection in a code such as NEKTON. 

We used Motif, which is the most popular, and, we believe, the most effective GUI 
system. Within this system we have adapted the current NEKTON user interface 
and added enhancements to achieve the above objectives. 

Many powerful tools are becoming available to allow the software developer quickly 
generate Motif GUI’s. Any description or list given here would no doubt be obsolete 
by the time this report reaches publication. 
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IMPLEMENTATION : 


The implementation of the new user interface features in NEKTON are done in the 
context of the Motif graphical user interface. Motif consists of 3 parts: the Mo- 
tif Window Manager, the Motif Interface Compiler, and the Motif Widget Library. 
NEKTON requires only the Motif Widget Library (libXm.a) and only when linking 
the graphical executables. Where this file is not part of the standard operating system, 
we ship precompiled and linked programs. Therefore, Motif requires no additional 
system resources for the new package. 

While the new code can and has been installed on VMS systems, in doing so we 
have concluded that new releases will not be supported on VMS. There are several 
reasons for this. First, the code relies heavily on the X Windows windowing system. 
This is, in theory, a subset of DecWindows and supported on VAXes. However, in 
practice, the performance of both the graphics and floating point operations on the 
CISC systems is low, and bog down the system. With the emergence of very powerful 
and low-cost RISC systems, the investment of time and money in VMS is dubious. 
For example, the DecStation 3100 is approximately lOx as fast as the MicroVax II, 
and the DecStation 5000 is in turn, almost a factor of 2 as fast the 3100. 

Moreover, while the VMS system has many advantages over UNIX, it is a different 
and proprietary system. Any enhancements taking advantage of these efforts can 
only be appreciated by a small group of users, and enhancements to UNIX (e.g., 
interprocess communication) would be curtailed to make it compatible with VMS. 
For all these reasons, it has been decided that VMS is an inappropriate platform on 
which a user should begin a committment to a new hardware/ software package. 

The use of Motif in the new interface has several benefits. First, interactive zoom and 
pan (scroll) are built-in to the geometry generator and postprocessor. The scrollbars 
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provide real-time “pulling” of the mesh, according to the standard conventions shared 
by Motif, Macintosh software, and other window- based programs. The “look and 
feel” of the program is customizable by the user via attributes files. 
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EXAMPLE CASES 


Mesh Generation 

The ability to quickly and easily generate meshes is important in MOCVD applica- 
tions. To facilitate this, we have incorporated the Motif pan/zoom capability into 
the PRENEK mesh generator. This enables quick and accurate input of meshes with 
elements with widely varying length scales, which are used, for example, in resolving 
thin boundary layers at the deposition surface of a rotating pedestal CVD reactor. 
In figures la and lb we plot the whole mesh and the mesh zoomed in on a corner at 
which we refine to get higher resolution near the corner. Thus, in entering a mesh 
we can at times focus in on specific areas, retaining the accuracy in inputting and 
modifying very small elements with the mouse. 

3-D Rotation 

In figure 2 we plot a three dimensional mesh for a 3-D CVD reactor undergoing 
rotation. This rotation is done via the slider bars, which allow for independent 
adjustment of pitch, roll, and yaw as defined by the user (screen) orientated frame of 
reference. This system allows for easy intuitive rotation of the mesh, as well as the 
full access to each of the 3 degrees of freedom of the mesh. Previous to this, only 2 
degrees of freedom were allowed. 

Contour Lines 

In figures 3a and 3b we plot the contours of temperature in a previous result of a 
CVD reactor undergoing g-jitter. This contour plot capability contains the same in- 
formation as the color fill plots, but is nonetheless useful for black-and-white displays 
and for submission of results in publications. 

Note also here that the zoom capability extends to all quantities and plot formats in 
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the post processor. 


Pulldown Menu Structure 

In figure 4 we plot the pulldown menu bar developed for POSTNEK. This menu 
bar structure enables easy access to all the capabilities of postnek. It minimizes the 
requirements for the user to read documentation on each function of the postprocessor. 
A scan through the POSTNEK menu bar and each of its submenus provides the user 
with an overview of the structure of POSTNEK and its capabilities. Moreover, this 
scan can be accomplished before committing to any action. In that sense, the process 
of selecting an item in a submenu is reversible. You can bring up the entire sequence 
of menus and submenus and decide at the bottom level whether to implement the 
action, or to abort the entire process without selecting the action. This is in contrast 
to more complex codes that present menus sequentially and only one level of the 
structure can be seen at a time. 

Equation Menu Structure 

In figure 5 we present the PRENEK set equation menu. This enables the various 
switches for the different aspects of the equation type to be switched on and off. 

Numerical Constant Input 

In figure 6 we show the Motif-style numerical input. This enables mouse-oriented 
editting of the default NEKTON constants. The standard click and double-click 
selection process, the cursor movement keys, and the delete operations work in a 
manner consistent with Macintosh and other mouse- oriented application programs. 

Scrollable Information Dialog Box 

In the upper left of figure 6 the dialog box logs informational messages from PRENEK 
and POSTNEK. The size of the vertical scrollbar in relation to its box is proportional 
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to the size of the visible text window to the total text. Thus, messages which have 
scrolled off screen are still available and can be retrieved by movement of the scrollbar. 
This enables large amounts of messages to be generated for the user without sacrificing 
excessive screen area. 

REFERENCES 

1. McMinds, Donald L., Mastering OSF/Motif Widgets, Addison Wesley, 
1992. 

2. The Open Software Foundation, The OSF/Motif Style Guide, Prentiss Hall, 
1990. 

3. The Open Software Foundation, The OSF/Motif User’s Guide, Prentiss 
Hall, 1990. 

4. The Open Software Foundation, The OSF/Motif Programmer’s Guide, 
Prentiss Hall, 1990. 

5. The Open Software Foundation, The OSF/Motif Programmer’s Refer- 
ence, Prentiss Hall, 1990. 

6. Nye, Adrian, Xlib Programming Manual Vll, O’Rielly and Associates, Inc. 
1990. 

7. Nye, Adrian, editor, Xlib Reference Manual Vll, O’Rielly and Associates, 
Inc. 1990. 


68 



1661 ZVS\ £ un f 


UIOJJ qS 3 J\ 1 U 3 UI 3 J 3 JBJJD 3 ds 


ixb :3urejsj uoissss 



■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■SBBSiiiiiSii 
■■■■■■■■■■■■■■■■■■■■■■■■ •■■■■■■■ ■•■■■■■■■■ 
■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■a 

■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■a 

ssssss8S88SES8:88S8S8888SSSsss8BS883ssess8sssssss8sb88sssss8 


■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■HHHJ 
■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■a 
'■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ ■■■■■■■■■■■■■■■■H.. 

■ ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■I 

iiiiiiimiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiRiiimiiiiiii 
■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■a 

.■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■a 

iiiaiiaiiviiiiiii>iiiiiiiiiiiiiiiiimiiiiiiiiiiiiiiiiiiiiiaiiaiiiiaiiiiiiiiiiiiiiiiRiiiiHiaiiH> 

■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■a 

BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB BBBBBiiBlBilgBIlBBBllBBglliiliBBBBBBBBBBBBBBBBSBBBBBaiBBiBlB 

IBBBBBBBBBBBBBBBB BBB BBB ■ BBBBBBBBBBBBBBBMaaaaaaaaMaaMaMaBMMMaaiiaaaaaBMMnHHHBMHHHBM 


IBBBBIBIBBBBB 

■■■■■»■■■■■■■ 

BBBBBBBBBBBBB 

BIBBBBBBBBBBB 


HHHHHbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbb 

BBBBBBBBBBBBBBBBBBBBBBBBBBBIB B1BBBBB111BB1B11 

M— 


BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB 

BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBg 

K ^HBBBBBBaiBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB 

bbbbbH 


■■BBBBBBBBBBBBB 

■■BBBBBBBBBBBBB 
■■■■■■■■■■BBBBBH 
■■■■■■■■■■■■■■■■■■■■I 


■■■■■■ ■■■■■■■■■■■■■■■■■■■■■■ 

■sssiiiiiissiaaaiiiMKSsssKSsaa^ 

iBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBgWBB BBBBBBBBj 

■■■■■■■■■■■■MMMMMMIBBMBRH 


■■■■■i 
PbbbbbM 

■■■■■■■■'& « m m m m m m m m w m m m m & m. m m l _ „ „ _ _ „ „ „„ _ .. 

Bbbbbbbbbbbbbbbbbbb| 
■MHHHHHHHHK B ■ ■ ■ BBB BBBBB B 
■HEBBBiiiiiBBBBBiiiBBBiBBBBBBBBBBBBBBBBBBM 
■■■■BBBBBBBBBBBBB BBB ■■■■■■■BBBBBBBBBBBBB ■BBBBBH 

■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■HI 

■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ BBBBB ■■■■■■■■■■! 


SSSSSSSSSSSSSSSSSSSSSSSESSSSSSSSSSSSSSSSSSSSSSBSSS.SBSBSSSSSSSSSSSSSSSSSS3SESSSSSSSSSSB8BB88BBBSBB 

SSHmHHHHHHHmmSm»SSBBBBBBBSBSBSBBBBBBB8BBSB8BSSSB8SBBBBSS8SBS| 

|j|SSSSBSBSSBBBSB8BSBS8BBSBBSS88SSSS8BSBSBS8BBSB 

■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■a 


■■■■■■■■■■■■BBBBB ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■^Hi 
■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■BBBBB 
■■■■■■■■■■■■■■■■■■■■ ■■■■■■ B ■■■ H ■■■■■■■■■■ B ■■■■■■■■■■ 

■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■a 
mmmmmmummmBmmmmumummummumummmmmmmmmammmmmmmmmmmmmmmm 
■■■■■■■■■■■■■BBBBB ■■■■■■■■! 


!■■■■■■■■■■■■■■■■■■■&■■ 
■HHHHHHHHKiiianBiiaiiiiiiBiiBBa 
■■■BBBaiSBBftaaBigBiBiBBi Sag ■!!■■■■■' ■■■■■■ ■■■■■■■■■■■ 
■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■HHMIHBH^H 

■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■I 

B ummm ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■] 
!■■■■■■■■■■■■■■■■■■■■■ ■■■■■■■■■■■! 
■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■I 
■■■■■ ■■■■■■■■■■■■■■■■■■■■■■■■■■■■I 

!■■■■■ BMMMHMMMHMMMMVHI 

■■■■■iiiiiiiiiiiiiiiiiiiaiiBBBiBil 

B ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■! 

■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■I 
[■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■a J 
■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■9 

■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■gg 
■■■■■■■■■■■■■■■■■■■■■■■ ■■■■■■■■B^l 
!■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ J| 
l■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■| 

■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■I 
■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■I 
!■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■! 
!■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■! 

[■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■I 
■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■| 

■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■I 
■■■■■■■■■■■■■■■■■■■■ ■■■■■■■■■■■■■! 
■■■■■■!■■■■■■■■■■■■■■■■■ ■ ■■■■■■■■] 


■■■■■■■■■■■■■■■■■■■■■■■■BBB 
IBaiBIIBBIIHIBIIBIIIBBBM| 
'■■■■■«■■■■■■■■■ ■■■■■■■ 
■■■■■■■■■■■■■■■■■■■■■I 

■■■■■■■■■■■■■■■■■■■■■HI 

!■■■■■■■■■■■■■■■■■■■■■■■■■■■ 
!■■■■■■■■■■■■■«■■■■■■■■■■■■■ 
■■■■■■■■■■■■■■■■■■■■■■■■■■■ 
■■■■■■■■■■■■■■ ■■■■■ ■■■■■■■■ 
■■■■■■■■■■■■■■ ■■■■■■■■■ ■■■■ 
■■■■■■■■■■■■■■■■■■■■■■■■■■a 
■■■■■■■ ■■■■■■■■■■■■■■■■■■■■ 
■■■■■■■■■■■■■■■■■■■■■■■■■■■ 
■■■■■■■■■■■■■■■■■■■■■■■■■■a 
!■■■■■■■■■■■■■■■■■■■■■■■■■■■ 
■■■■■■■■■■■■■■■■■■■■■■■■BBB 
■ ■■■■■■■■■BBMBBBBBBBIBBBBB 
>■■■■■■■■■■■■■■■■■■■■■■■■■■■ 
■■■■■■■■■■■■■■■■■■■■■■■■■■■ 

{■■■■■■■■■■■■■■■■■■■■■■■■■■■HBI 

mmmaMmmmmmmaammttammaamaammrnw 


sTnnnnnnnssHHHHSHHHHHSHSHH 

UjySj|8S8888S888SS888S8S885s8888SSS88SS8SSSS 

Isess:::::::s88se:s::s:sss8888SS88S8S88 

8S8SSSE88SS88S88888888SS888S8S888S8S8SS8SS8SS 

■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■SB 
■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■I 

■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■a 
■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■! 
■■■■■■■■■■ BBBBB ■■■BBBBB ■■■■■■■■■■■■■■■■■£■■■■ 


■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■a 

■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■a 

■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 


!■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 

fiaasa8as85SSH5aaasaaHsssssssssssssssassssi 

I888SS8SSE8SSS88S8SS8! 
■■■■■■■■■■■■■■■■■■■■< 
■■■■■■■■■■■■■■■■■■■■I 

■■MM ■■■■■■ ■■■■■■■! 

■■■■■■■■■■■■■I 

iiiiiii«B ■■■■■■ ■■■■■! 

■■■■■■■■■■■■■■■■■ilij 

RlH ■■■■■■■■■ 

■■■■■■■■■■■■■■■■■ 
■■■■■■■■■■■■■■■■■ 
■■■■■■■■■■■■■■■■■iiS 

ssssssssSssIssisiSni 
■jjjjjjs888S888SSSjj 

^^^VHbbbb ■■■■■■■■■ 

■ ■■■■■■■■■■■■■■■■■■■! 

■■■■■■■■■■■■■■■■■■■■j 

■ at ■ ■■■■■■■■■■■■■■■■■! 
■■«■■■■«■■■■■»■■■■■■ 
■■■■■■■■■■■■■■■■■■■|; 
■■■■■■■■■«■■■■■■■■■■ 

:::::::::::::::::::: 

SSSSSSSS^ ■■■■■■■■■■! 

■ ■■■■■■■■■■■■IIIiHH 





Figure 2 






NEKTON 2 : 



72 








73 


Session Name: 
Postnek Results 
Jun 3 14:48 1991 



Figure 4 




















teem Vartlan 2.7 
Qigggi i Nmw fur This Sssslon: 
j^lwuw S— >lon prcbl«rf 
mil had PrmUra Froo praklMl.r« 

01800 SETTING: STEMV TMD ■INENSIOWL WftT TTWf 
OMBO SETTING: IMSIEMT HO ■IIWSIWflL HEAT HI 

cmot SEniNG: unstewy no micmsiqm. sums 
atm setting: unstemt no bihekioml nrwiw- 

Tim ftwplnfl scHm* mC to Flrot (Ww 

amnsniNC: unstemy no bdosioml nwioi- 


f l 


ii 


MMM27 













76 



























ZHIL . 


Chapter 6 


Low Mach Number 
Compressibility 


OBJECTIVE: 

The objective in this report period was to implement the low Mach number compress- 
ibility into NEKTON and demonstrate this capability for CVD problems. We have 
shown the use of this accurately modelling the effects of compressibility due to the the 
effects of varying temperature. This is an improvement over the Boussinesq approx- 
imation which was previously used to model the effects of density changes. In this 
report period we have implemented this capabililty in NEKTON for general 2-D, 3-D 
and axisymmetric geometries. We have demonstrated this new capability in simple 
geometries with which we have verified our results through comparison with analytic 
solutions. We have also demonstrated this capability on realistic CVD geometries. 

BACKGROUND: 

Solution of the full time-dependent compressible Navier-Stokes equations presents 
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difficulties in combining the vastly differing acoustic and viscous time scales. In this 
approach we are able to decouple the acoustic compressibility, with its numerical and 
physical complexity which is irrelevant to CVD applications. We retain the thermal 
compressibility term which is crucial to accurate solution of CVD problems. 

CVD problems typically involve absolute temperatures which can vary by a factor 
of two. The Mach numbers involved in CVD are typically of order 10 -3 to 10 -4 . 
This combination of parameters verifies that there is a need for a scheme which will 
account for the density variations due to temperature, and need not account for 
density changes due to momentum. 

Torczynski [3] studied the inviscid problem in which a perfect gas is heated by means 
of energetic particles. He was able to confirm the validity of the acoustically filtered 
equations of motion for low Mach number by comparison with solution of the full 
gasdynamic problem. However, the inviscid assumption is valid when used as in most 
high Re aerospace applications where viscosity affects only thin boundary layers near 
walls or thicknesses of shock waves. 

CVD applications, however, involve significant viscous and diffusive effects, and there- 
fore require solution of the acoustically-filtered full Navier-Stokes equations. Addi- 
tionally, as Makarov [4] noted, it is crucial to include the additional viscous and 
diffusive term associated with thermal gas expansion (“bulk viscosity”) that arises 
in the compressible case. This is because in most CVD geometries a boundary layer 
forms normal to the deposition surface, and this viscous compressible term is signifi- 
cant in quantifying this boundary layer thickness and the resultant deposition rates. 

This formulation offers significant advantages over solution of the full compressible 
Navier-Stokes equations. It greatly reduces the computational load presented by the 
presence of the viscous and convective time scales of order seconds and the acoustic 
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time scales of order fractions of milliseconds. With this formulation it is possible to 
concentrate on study of the real CVD problem without introducing the complexity 
of irrelevant and short-lived acoustic waves. 


COMPRESSIBLE TRANSPORT EQUATIONS: 

The development of compressible transport equations for low Mach numbers 
(characteristic of CVD systems) stems from the need to fully capture the effect 
of temperature gradients in the system on the flow field m the reactor. Subsonic 
compressible flow models have been reported in [1,2, 4-7]. More recently Strelets and 
co-workers [8] have extended this formulation to concentration convection and multi- 
component chemically reactive flows. The subsonic model takes advantage of the low 
Mach number of the system to eliminate acoustic waves and results in a system of 
equations similar to that of incompressible transport equations. 

Non-dimensionalizing the length, velocity , time, pressure, temperature, and 
properties of the gas: 


x,y,z, V,t,P,T, p,fi,k 

by system reference values 

L\ V\ L*/V*, p*RT *, r , />*, p\ **, 

respectively, results in the following equations: 

0 


( 1 ) 
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(2) 


DV 
' Dt 


1 V7 n.sL* . , 1 

^p VP+ V* pe ° + Te T 


> 1 * = + ^ V -( Wr ) + i(^ l)(M 2 )[2MW + (VK)‘)^V.^I] (3) 


P — pT (4) 

t = ,4W+(W) t -Hv.KI] (5) 

o 

In the above g is the gravity, 7 the ratio of specific heats of the gas, M the 
Mach number = V*/c*, e g unit gravity vector, Re the Reynolds number = p*V*L m /p m , 
Pr the Prandtl number, and I the identity matrix. Also the superscript T denotes 
the transpose of the operator. 

The subsonic model is based on two assumptions which are valid in CVD 

systems: 

1. The Mach number M = V*/c* << 1. 

2. The hydrostatic compressibility e = gL*/RT* « 1 

In order to derive the subsonic compressible equations the field variables are 
expanded in terms of 7 A/ 2 and t. For example, pressure is expanded as: 


P = p(°) + 7 m 7 pw + € pm + .... 


( 6 ) 


80 



Inserting these expansions to eqs. (l)-(5), equating like-powers of the two perturba- 
tion parameters, and dropping the (0) superscript for all variables except pressure 

one obtains: 


V P (0) = 0 

(7) 

% + P V.V = 0 

(8) 

= -(vp ( 10) + + f £pe, + jb(r) 

(9) 

p EL _ 7 “ 1 dp{0) + 1 v-(jfcvr) 

p Dt 7 dt RePr 

(10) 

pW = P T 

(11) 


By recognizing that in equilibrium V P< 01 > = p.e a , where p , is the equilibrium density 
of the gas, eq. (9) can be rewritten as: 

P^ = -Vi ><10) + ^(p - P.Y, + J^( T ) < 12 > 

Equations (5)-(7), and (9)-(ll) imply that in the subsonic model the pressure 
term splits into two terms: a thermodynamic pressure P {0) which is only a function 
of time (from eq. (7)) and enters the equation of state, and P (10) which is function of 
both space and time and accounts for the dynamic effects. In the subsonic formulation 
the viscous heating term of the energy equation is also dropped. By differentiating 
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the equation of state and using the above equations a new form of the continuity 
equation 'is obtained: 


1 1 1 dPW 

< 13 > 

For open systems the pressure term P^(t) is constant and one obtains the 
following set of equations: 



(14) 

DT 1 „ ,, 


>l»=RePr V -^ 

(15) 

VV = ^RePr V -^ 1 

(16) 

P (0) = pT 

(17) 


The above sets of equations (for open CVD systems) is very similar to that of 
incompressible transport equations which are currently solved (in three space dimen- 
sions) by NEKTON. NEKTON hats been modified to solve the above set of equations 
in order to correctly simulate the effect of local thermal expansion of gas during 
CVD processes. This reformulation of the governing equations takes advantage of the 
spectral-element discretization and iterative solvers already implemented in NEK- 
TON. 


82 



IMPLEMENTATION : 


The above procedure was implemented in NEKTON in the context of the 
current formulation for solving the elliptic equations associated with incompressible 
flow. The current implementation involves solving two problems: The Navier-Stokes 
equatons for viscous incompressible flow, and the energy equation for the temperature. 
The coupling of the convective terms was achieved by iterating between these two 
calculations in the process of time-stepping. 

An essential difference in the calculation of low Mach number compressible 
flow in the enforcement of the divergence condition of the velocity. In the Navier- 
Stokes calculation the incompressibility condition enforced with V • V = 0 is replaced 
with a condition (V • V = -CDT/dt) relating the divergence to the total derivative 
of the temperature field. In this manner the methodology of the solving the (ellip- 
tic) system of equations guarantees that the pressure is solved globally in a coupled 
manner, without any local mechanisms which cause propagation of pressure (sound 
waves). 

Additionally, modifications were made to the viscous stress tensor to include 
the extra viscous divergence terms and finally a body force proportional to density 
was added to replace the Bousinnesq approximation. The significant coding tasks 
involved in these changes were to implement and verify the above equations, and 
extensive changes from an assumed constant density to a variable density throughout 
the code, including the matrix preconditioners, etc. 

We were able to introduce this capability in the context of the full code 
NEKTON. Thus this new capability can be used in two-and three-dimensional as 

i .~Z . : - ... 

well as axisymmetric geometries. Additionally, this compressibility can interact with 
the large set of boundary conditions available in NEKTON. 
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EXAMPLE CASES 


In this section we present two classes of examples. In the first we present 
simplified geometries amenable to analytical analysis in order to demonstrate the 
validity and accuracy of the implementation. In the second class of problems we 
demonstrate the use of this method on CVD applications and compare the results of 
the compressible solution with that calculated using the Boussinesq approximation. 

TEST GEOMETRIES 

In the first we analyze the one-dimensional unsteady problem described in figure 1. 
Here we begin with a flow whose velocity, temperature, and pressure are unity. We 
set the Mach number to be 0.1, the Reynolds number to 100, and the Prandtl number 
to 1. We fix the inflow velocity and temperature to be unity at the left end of the 
domain and set the right edge of the domain to be outflow conditions with a pressure 
perturbation fixed at 0. Note that all pressures plotted are the perturbed (Pi) values; 
the constant pressure P 0 is fixed as input. Symmetry conditions ( d/dn = 0) are used 
at the top and bottom boundaries. The thermal conductivity k = 10 -3 , viscosity 
v = 10 -3 , and domain length of 1. We input a heat flux of q=2k + 2/R (1-x). We 
expect an exact solution at steady state of 


T = 1 + 2x — x 2 
V = 1 + 2x — x 2 
p — 8/3//(l — x) = 2.66xl0 _3 (l — x) 

In figures 2-4 we plot the steady values of velocity, pressure, and temperature, respec- 
tively. In figures 5-7 we plot the time evolution of these quantities at the midpoint 
of the computational domain. In figures 8-10 we plot quantitatively the profiles of 
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these quantities along a line from inlet to outlet. At the midpoint of this line in 
figure 9 a small change in slope of the solution appears where the two elements of the 
computational domain meet. This error is associated with the finite spatial resolution 
and can be made arbitrarily small by increasing that resolution. In figure 11 we plot 
the Li error in these quantities as a function of time step magintude. These results 
confirm the accuracy of the solution and confirm that the scaling of the error with 
time step is linear, as expected in our first order stepping scheme. 

In figure 12 we plot the geometry of a two-dimensional channel flow past a heated 
section of otherwise adiabatic walls for a Mach number of 0.1. In figure 13 we plot the 
streamlines of the steady flow past the heated section. In figure 14 we plot contours 
of the speed. Note the increase in speed as the density decreases downstream of the 
heat source. In figure 15 we plot perturbed pressure (Pi) contours. Note that the 
magnitude of this scales correctly (as Ma 2 ) confirming a one percent expected error 
(Pi/P 0 ) for Ma=0.1. In figure 16 we plot the temperature contours. In figure 16b 
we plot the coutours of divergence. Note the peaks in the divergence of steady flow 
are near peaks of temperature gradient, confirming the expected behavior of the flow 
divergence scaling with DT/dt. Finally in Figure 17 we plot the temperature history 
downstream of the heated plate, which shows the approach to steady state. 

CVD REACTOR GEOMETRY 

In figure 18 we plot the geometry of an axisymmetric CVD reactor. We use the 
following properties of Hi and flow parameters: 
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#2 Properties 


R 

C P 

7 

Simulation Parameters 


Q 

= 

10 Std L/min 

Po 

= 

0.1 Atm 


= 

1.01 x 10^ 

m i 

T xn 

= 

300 °K 

*u>a fer 

= 

1000° A' 

T 

= 

300 + 1000 = 650°* 

c 

= 

\JlRT = 1945 m/sec 

D 

= 

0.1m 

Vinlet 

= 

0.2m /sec 

M 

= 

rH 

o 

o 

o 

o 

II 

Re 

= 

pVD _ 0.37 x .2 x . 
p ~ 1.5 x 10- 5 

7 

= 

//f = .0037fc<?/m 3 

Pr 

= 

0.7 

k 

— 

= .307 " 
Pr m°K 


= 1.5x10 


-5 


Nt-s 


m‘ 


= 8134 
= 4157 


= 14315 


= 1.4 


kg — mol —° I\ 

J 

kg K 

J 


1 kg — mol 

2kg 


kg I< 


( 6 . 1 ) 


( 6 . 2 ) 
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Note that the low Mach number of 10~ 4 justifies our assumption for the accuracy of 
the low Mach number assumption for CVD applications. We compare cases using the 
incompressible code and that with the low Mach number compressibility. In figures 19 
and 20 we plot the velocity vectors and temperature contours using the incompressible 
assumption with properties based on the average temperature. In figures 21 and 22 
we plot these same quantities using the low Mach number approximation. Note the 
large difference in peak velocity (1.2 vs. 0.5), resulting from the large volumetric 
expansion of the gas from 300 K to near 1000 K. 

CONCLUSION 

We have demonstrated that the low Mach number approximation is appropriate for 
CVD applications. In this report period we have concluded work on the develop- 
ment and implementation of the low Mach number scheme in NEKTON. We have 
demonstrated this new capability on theoretical problems and on problems in CVD. 
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Chapter 7 


Point-Implicit Integration 


OBJECTIVE: 

The objective in this report period was to implement an implicit time integration of 
the chemical reactions into NEKTON. This was to be done using the point-implicit 
method described by Bussing and Murman [1]. 

In this report period we have implemented this capabililty in NEKTON for general 
2-D, 3-D and axisymmetric geometries. We have demonstrated this new capability 
in simple geometries with which we have verified our results through comparison 
with analytic solutions. We have also demonstrated this capability on realistic CVD 
geometries. 

BACKGROUND: 

In many cases, solution of chemically reacting flows is made difficult by wide sepa- 
ration in the time scales involving the diffusive and convective fluid effects and the 
rates of chemical reactions. In typical MOCVD problems, the reactions can occur in 
time scales which are orders of magnitudes faster than the fluid effects. The physical 
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interpretation of this is that the reactions are diffusion-limited, that is, the reactants 
are quickly depleted locally and the reaction is governed by the slower time scales 
associated with transport of new reactants to the reaction site. These reactions re- 
sult in a stiff numerical problem, with the stiffness usually defined by the ratio of the 
smallest time scale to the largest time scale, 


StlffnCSS — ^largest! ^"smallest 

Bussing and Murman [1] developed a point-implicit method for use in solving these 
stiff problems for use with Euler equation calculations using finite volume methods. 
We have adapted this point- implicit method for use in solving the compressible and 
incompressible Navier- Stokes equations in NEKTON. 

EXAMPLE CASES 

To test the consistency, time-accuracy, and accuracy in the steady-state solution, two 
test problems were chosen. These were chosen to be relatively simple problems which 
had analytic solutions, but were nonetheless representative of calcualtions required 
for MOCVD. Each reaction has the form of: 


R = kS 


where the reaction rate constant k for a species of concentration S was given values 
of 1 and 100. This problem was run in a reaction/diffusion environment described 
by: 


DS 

Dt 


-kS + V 2 S 
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For this the steady solution is of the form 


S = Cie 10x + C 2 e~ l0x 


and 


S = C,t x + C 2 t~ x 


TEST RESULTS 

These calculations were run in the domain x=[0,l] with S(0,t)=l; S(l,t)=2; and S(x,0) 
= 1.0. the resultant steady solution is: 

S = 0.6944e r + 0.3056e- e 
5 = 9.07978E - 5e r + 0.999909202e _x 

The steady concentration solution calculated by NEKTON corresponding to the fast 
and slow reactions is plotted in Figures 1 and 2, respectively. 

The error in the steady result (as measured by concentration at the domain midpoint) 


is summarized as follows: 


Time Step 

Fast Reaction 

Slow Reaction 

SS Soln 

SS Error 

SS Soln 

SS Error 

0.001 

1.33023 

0 

0.0202140 

0 

0.01 

1.33023 

0 

0.0202140 

0 

0.1 

1.33024 

10- 5 

0.0202128 

10- 6 

1.0 

1.33024 

10" 5 

0.0202128 

10- 6 

10.0 

1.33024 

10- 5 

0.0202122 

10- 6 

100.0 

1.33024 

10" 5 

0.0202122 

10* 6 
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The time histories of the approach to the steady solution is plotted in figures 3a-4c. 
Figures 3 a,b, and c represent the history of the concentration in the fast reaction 
as it approches its steady value with time steps of 0.001, 0.01, and 0.1, respectively. 
Figures 4a, b, and c represent the same for the slow reaction. 

DISCUSSION 

In as much as it is possible, this scheme preserves the time accuracy of the algorithm. 
For those component reactions whose characteristic time r is small compared to the 
time step, a time-accurate calculation of its contribution is done. For those parts 
of the simulation whose characteristic time r is large compared to the time step, a 
stable iteration occurs which leads it to the correct steady state. This process results 
in a coupled calculation in which the faster processes are solved in a quasi-steady 
iterative scheme, and the unsteadiness is governed by the time-accurate calculation 
of the slower system components to which the fast process is coupled. For problems 
in which the time step is much greater than the characteristic time (e.g., DT=0.1 for 
a reaction with time scales of 0.01, figure 3c), the time stepping algorithm effectively 
becomes an iteration scheme, with the species converging to steady state in a few 
time steps. The advantage is in that many less iterations are required in this pseudo 
time-stepping scheme than actual time steps required in the standard time-stepping 
scheme. 

CONCLUSION 

We have demonstrated that the point-implicit method in NEKTON for solution of 
stiff chemically reacting flows. It has been demonstrated to give full time-accurate 
solutions for those problems which are temporally resolved, that is, for problems in 
which the time step is small enough to capture the time scales for the physics of 
interest. It has been demonstrated to reach accurate steady state solutions for those 
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problems that are not temporally resolved. Where both fast and slow time scales 
exist in a coupled problem, the enhanced NEKTON gives time-accurate solutions for 
the slow components and smears the fast reactions across a few time steps. 
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Non-dilute Mixtures 
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OBJECTIVE: 


The result achieved in this report period was to implement the effects of non-dilute 
mixtures in the convection-diffusion equations in NEKTON. 

BACKGROUND: 

The chemistry in the gas phase in CVD reactors and the concentration of each 
species of mass involved in the reaction is important in finding the overall rates of 
deposition and the spatial uniformity of that deposition. The standard diffusion equa- 
tions used when small concentrations of contaminants diffuse through background 
media passively (without affecting the background media) break down when the con- 
centrations become large. 

IMPLEMENTATION : 

We begin with the relations for multicomponent nondilute mixtures developed in 
White (1). We have a mixture of n fluids, with individual total masses m 3 , ... 

contained In volume v. The density of this mixture would be defined as: 


where 

and 


m 

P = — 


m = 2 m, 
1=1 


rrti 

Pi = — 

v 


(i) 


( 2 ) 


( 3 ) 
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where 


P = Y,Pi 

i=i 


we define u as the mass fraction: 


which sum to one: 


The velocity V is defined as 


m p 


- 1 

t=i 


( 4 ) 


( 5 ) 


(6) 


V = -X>V t = £u;,V, 
P i={ i=i 


( 7 ) 


where p{V \ = f the mass flux of species i. 

The difference between the mass averaged velocity and the species velocity of 
component i is the diffusion velocity of component i, (V; — V). The diffusion mass 
flux is: 


j,- = MV.- - v) 


( 8 ) 


The net diffusion mass flux sums to zero: 

= ° 


( 9 ) 


1=1 


The implications of this in implementation of these equations in NEKTON is as 
follows. The boundary condition for absorption of a given species is p,V, = f,. The 
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boundary condition for a given species that is not absorbed is piVi = f, — 0. In the 
vector quantities used here and in what follows, the coordinate system used is local, 
with all vector quantities (fluxes, velocities) represent components oriented normal 
to the surface with the convention that quantities oriented outward from the domain 
are positive. From (7), the velocity at the surface can be calculated as: 

v = -£f, = -X>v, ( 10 ) 

pf^i p i=i 

Note that in cases in which one or more species are being absorbed, i.e., f, ^ 0, 
the velocity at the wall is nonzero. 

The mass balance and associated boundary conditions for the transport of chem- 
ical species is, then: 

Knowing V and V,', The diffusion mass flux at the surface is. 


j, = p, (Vi - V) = D,Vp, 


(ID 


This mass flux can be input as a surface flux. 

Continuity Equations: 

The continuity equation for a single component gas under low Mach number com- 
pressibility is 


1 dT 

V • V = — — 

T dt 


( 12 ) 


With a mixture of ideal gases the effective gas law constant for the mixture is: 

7?=i (13) 

pU 
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The continuity equation for the gas mixture under low Mach number compress- 
ibility is 


1 dT 1 dR _ 1 dT 1 A n d Pi 
T dt + H dt T <ft + PiRih'i dt 


( 14 ) 


while the momentum equation becomes: 


DV 

Dt 


VP 

— + vV 2 \ + / 


( 15 ) 


SCALINGS: 
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For a binary mixture of gases a and b in which a is being absorbed and the b is 
not, we have the following: 


V = 




The diffusion mass flux of b at the surface is: 



j* - A(V* - V) = D b V Pb (17) 

Since Vj, = 0, we have 


-PbV = D b Vp b (18) 

_ T * ! {t t •. * ...• . ; 

Fo a boundary layer thichness 8 , the conentration gradient Vp b scales as 1/8, so 
that 

(19) 

_ v_ 
— D 
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e.g., the dimensionless boundary layer thickness^d/Tjscales as where Sc 
and ReSc — 



Chapter 9 


Project Summary 


CHEMICAL VAPOR DEPOSITION 
FLUID FLOW SIMULATION 
MODELLING TOOL 

Chemical Vapor Deposition (CVD) processes combine problems from several 
disciplines. They require solution of the flow problem to model the transport of heat 
and chemical species, as well as solution of the chemical kinetics problem within the 
flow. Modeling of the of the flow field requires solution of the Navier-Stokes and energy 
partial differential equations which are coupled through thermal expansion, buoyancy, 
and convection. The principal effect of the flow field on deposition is that it affects 
the transport of the species to the reaction sites and their temperatures. Modelling of 
the chemistry consists of solution of a set of ordinary differential equations governing 
the (local) reaction chemistry. 

To these sets of problems, we have applied the spectral element method 
(SEM) using the code NEKTON. The SEM combines the methodology of the standard 
finite element with high-order (spectral) techniques. In this way the accuracy of 
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spectral techniques can be combined with the geometric flexibility of finite element 
methods. 

NEKTON incorporates additional innovative features in its design that main- 
tain it at the leading edge of the state of the art in numerics. It is the first commer- 
cial CFD code to be intrinsically designed for use on parallel processors. It is, to our 
knowledge, the only major commercial CFD package which can effectively utilize par- 
allel processors. Its locally structured, globally unstructured spectral element mesh 
is naturally suited to a geometry-based parallelism in which each spectral element 
(or group of spectral elements) is mapped to a separate processor/memory, with the 
individual processor/memory units linked by a relatively sparse communications net- 
work. This structure makes for high parallel efficiency. Another innovation is that it 
is designed to be inherently three-dimensional, so that features implemented in two 
dimensions are automatically available in the three-dimensional code. 

In modifying NEKTON to address issues in CVD modelling, we have made 
progress in three essential areas: First, e have implemented a number of features 
into the code which aid in the analysis of CVD problems. Second, we have devel- 
oped improved and robust numerical methods. Third, we have identified methods of 
programming methodology that will result in future codes that are more portable, 
modular, and maintainable. 

The capabilities implemented were specific to CVD and especially to CVD in 
microgravity applications. We first implemented the capability of solution of problems 
in which the “g-jitter” typically found in microgravity applications was a significant 
influence on the convective transport. This was done through user specification of an 
arbitrary time-dependent and orientation dependent gravity vector. 

The major task in this project was the implementation of low-Mach num- 
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compressibility. Whereas most analyses of CVD processes assume a Boussinesq ap- 
proximation for thermal expansion, Makarov and Zhmakin [1] have recently shown 
that this approximation can result in errors of the order of 100%. In contrast, we 
have shown that the low-Mach number approximation gives errors of less than 1% 
through the range of CVD problems. This approximation takes advantage of the fact 
that the compressibility in CVD applications is essentially due to thermal expansion, 
and that compressibility due to the momentum equation is negligible. This gives the 
important advantage in solution of CVD problems in reduced physical and numeri- 
cal complexity. Additionally, the stiff numerical problems associated with the widely 
disparate time scales of the fully compressible Navier-Stokes equations are avoided. 

Another important feature for CVD was the ability to solve problems in 
which the properties (thermal and mass diffusivities) vary with temperature and 
species concentration. 

A feature crucial for rotating pedestal CVD reactors is the ability to include 
swirl (azimuthal velocity) in the calculations. The swirl calculation enables rotational 
pumping of the flow by the rotating pedestal to be included in the model. The swirl 
velocity is coupled to the planar axisymtric velocities in the reactor via the coriolis 
and centrifugal forces. 

An additional feature has been implemented into the code. While the code is 
inherently iterative and unsteady, we have identified some CVD problems for which 
a steady, direct solver is more appropriate. 

We have also tested modules for multi-body radiation, Soret effects, and 
nondilute mixtures. We have to date included into the commercial code several 
numerical capabilities developed in this project. The major developmental goal of 
this project, compressible transport in the gas phase, has been fully integrated into 
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the code. We have also incorporated single-body radiation, the steady solver, the 
ability to put in arbitrary and time- and orientation-dependent gravitational forces 

in order to model g-jitter. 

The programming advances that resulted from this project were as follows. 
The implementation of the direct solver presented memory requirements that were 
substantially larger than the iterative solvers. Moreover, the memory for the direct 
solver varied as square of the the number of degrees-of-freedom, rather than linearly, 
resulting in huge variations in memory requirements for different problems. We have 
developed the capability in our commercial code of mixing C routines with the FOR- 
TRAN which makes up the bulk of the code. The C code functions to allocate memory 
during run-time to pass to the FORTRAN subroutines. This approach was so suc- 
cessful and attractive that we are currently in the process of a complete rewrite of 
the code into the C language. Additional benefits of C are modularity, encapsulation 
of data, and natural transition to the object-oriented language C++. 

In the Phase III commercialization of this code we expect to incorporate the 
remainder of the features developed unde r this pr og ram i nto the commercial NEK - 
TON. Using the newly developed progrTmming tecTinlqu^ We fexpect the maSieS* 
modularity of the code to enable us to re-use modules and retain a broad range of 
features in a rapidly evolving code. 
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